Outline

• Spline approach
  • Natural cubic spline
  • Regression spline
  • Smoothing spline
  • Example: B-spline
    
  • Example: Principal Component Analysis
    
• Local polynomial approach
• Curse of dimensionality

Natural Cubic Spline

• How to find a good-looking curve to interpolate all points $(x_i, y_i)$ ?

Mathematician says… find a curve that minimize its second derivative over whole range of input x

$$J(f)=\int_{x_1}^{x_m} f''(x)^2 \text{d}x$$

• Theorem

Of all continuous function that are continuous on $[x_1, x_m]$, have absolutely continuous first derivatives and interpolate $\{x_i, y_i\}$ natural cubic spline $g$ is the one that is smoothest in thes sense of minimizing $J(f)$

Natural Cubic Spline

• Definition (natural cubic spline)

Consider $\{x_i,y_i:i=1,\dots,m\}$ where $x_i<x_{i+1}$ the natural cubic spline $g(x)$ satisfies

1. on each interval $[x_i,x_{i+1}]$, $g$ is a function made up of sections of cubic polynomial.
2. $g$ interpolate points $\{x_i, y_i\}$. i.e., $g(x_i)=y_i$
3. except for boundary points, the whole spline is continuous to second derivative.
4. $g''(x_1)=g''(x_m)=0$

• an important idea: With $k$ knots ($k=m$), we use $k-1$ peices of cubic polynomial to construct a spline function.

Regression spline

Consider defining a cubic spline function $f(x)$ with $k$ knots, $x_1,\dots,x_k$.

Let $\beta_j=f(x_j)$ and $\delta_j=f''(x_j)$. Then

$$f(x) = a^-_j(x)\beta_j + a^+_j(x)\beta_{j+1} + c^-_j(x)\delta_j + c_j^+(x)\delta_{j+1} ~~~ \text{ if } x_j\leq x\leq x_{j+1}$$ where $$\begin{align*} a_j^-(x)=&(x_{j+1}-x)/h_j \\ a_j^+(x)=&(x-x_j)/h_j \\ c_j^-=&[(x_{j+1}-x)^3/h_j-h_j(x_{j+1}-x)]/6\\ c_j^+=&[(x-x_j)^3/h_j-h_j(x-x_j)]/6 \\ h_j=&x_{j+1}-x_j \end{align*}$$

• I know, it is disgusting. Fortunately, there is anoher kind and beautiful representation of spline function $f$ $$f(x) = \sum_{i=1}^k b_i(x)\beta_i$$

Regression spline

$$f(x) = \sum_{i=1}^k b_i(x)\beta_i$$

• This formulation is important and frequently to be used in the research, we call it basis expansion
• In basis expansion, it’s not necessary to set number of knots ($k$) to be equal to number of observation ($n$). $k$ could be smaller than $n$. Note that when $k<n$, the fitted smooth curve will not interpolate all data points.
• $b_i(x)$ is a prespecified basis, till know, we use cubic spline basis only.
• Any kind of basis that is rich enough can be used to represent the true function $f$. It is similar to build up the same castle with different Lego bricktopias.

• Other rich basis
  • Truncated power basis
  • B-spline basis
  • Wavelet Haar basis
  • eigen basis (Principal component analysis)

Smoothing spline

Statistican says that we just need an approximation and in addition, we’d like to control the roughness

$$\sum_{i=1}^n\{y_i-f(x_i)\}^2 + \lambda\int f''(x)^2\text{d}x$$ where $f(x)=\sum_{j=1}^kb_j(x)\beta_j$

• What kind of basis should we use ?
  • It’s a long story… different basis has its specialty, it depends on you purpose.
    
  • Basically, we can see that cubic spline basis overlapes with its neighbors. If the overlapes is servere, the coefficient $\beta$ may be difficult to be estimated.
    
• The choice of tuning parameter $\lambda$ is an important issue, it can be tuned by CV, GCV or AIC… etc.
• Compared to the choice of $\lambda$, choice of knots is not that important. In the literature, people suggest that we can use many knots to do basis expansion and control $\lambda$ to adjust the roughness of penalty.

Example: Smoothing spline with B-spline

$$\sum_{i=1}^n\{y_i-f(x_i)\}^2 + \lambda\int f''(x)^2\text{d}x$$ where $$f(x)=\sum_{i=1}^k B_i^m(x)\beta_i$$ $$B_i^m(x) = \frac{x-x_i}{x_{i+m+1}-x_i}B_i^{m-1}(x) + \frac{x_{i+m+2}-x}{x_{i+m+2}-x_{i+1}}B_{i+1}^{m-1} ~~~ i=1,\dots,k$$ $$B_i^{-1}=I\{x_i\leq x\leq x_{i+1}\}$$

library(mgcv)
set.seed(2) 
n <- 400
dat <- gamSim(1,n=n,dist="normal",scale=2)

Gu & Wahba 4 term additive model

x=dat$x2; y=dat$y
plot(x, y)

###
### B-spline
###
library(fda)
bsp <- create.bspline.basis(norder=3, breaks=place.knots(dat$x2, 39))
B <- getbasismatrix(dat$x2[order(dat$x2)], bsp, returnMatrix=TRUE)
D <- diff(diag(40), diff=2)
lambda <- 1
a = solve(t(B)%*%B + lambda*t(D)%*%D)%*%t(B)%*%y[order(x)]

plot(x, y)
for(i in 1:ncol(B)) lines(x[order(x)], a[i]*B[,i], type="l")

plot(x, y)
for(i in 1:ncol(B)) lines(x[order(x)], a[i]*B[,i], type="l")
lines(x[order(x)], B%*%a, lwd=2, type="l")

lambda <- 30
a = solve(t(B)%*%B + lambda*t(D)%*%D)%*%t(B)%*%y[order(x)]
plot(x, y)
for(i in 1:ncol(B)) lines(x[order(x)], a[i]*B[,i], type="l")
lines(x[order(x)], B%*%a, lwd=2, type="l")

Example: principal component analysis

Do you remember the Taylor expansion ? expand $f(x)$ at point $x=x_0$

$$f(x) = f(x_0) + f^{(1)}(x_0)(x-x_0) + \frac{f^{(2)}(x_0)(x-x_0)}{2!}(x-x_0)^2 + \cdots$$

The former term describes the curve near $x_0$ (more important), the latter terms away from $x_0$ (less important).

Example: principal component analysis

Actually, there is a functional version… Karhunen Loeve Decomposition

Consider a centered stochastic process ${X_t}, t\in[a, b]$ satisfying a technical continuity condition and with covariance function $K_X(s,t)$. $X_t$ admits a decomposition

$$X_t = \sum_{k=1}^{\infty} Z_ke_k(t)$$ where

• $e_k$ be an orthonormal basis of $L^2([a,b])$ with respect to a eigenvalue $\lambda_k$
• $Z_k=\int_{[a,b]}X_te_k(t)\text{d}t$
• $E[Z_k]=0 ~\forall~ k\in\mathbb{N}$ and $\text{var}(Z_k)=\lambda_k$
• $\sum_{k=1}^{m} Z_ke_k(t)$ converges to $X_t$ in $L^2$ and uniformly in $t$ when $m\rightarrow\infty$

• It can be shown that KL is the best orthogonal basis upon the process in the sense of minimizing $L^2$ norm

• Is it familiar for you with the basis expansion?

  $$f(x) = \sum_{j=1}^k b_j(x)\beta_j$$

Smoothing in high dimension

• Thin plate spline basis

• Tensor product approach

Example: image reconstruction (Eigen basis)

• In this example, we use tensor product approach with eigenbasis to reconstruct an $658 \times 718$ image.

library(biOps)

std <- function(x){
  y <- matrix(c(0,255),2,1)
  A <- cbind(range(x),1)
  coefs <- solve(A,y)
  xmat <- cbind(x,1)
  z <- xmat%*%coefs
  z <- round(z)
  z[z<0] <- 0
  z[z>255] <- 255
  z
}

y = imgGreenBand(readJpeg("Ansh.jpg")) # 658 718

s = svd(y)
u <- s$u
v <- s$v
d <- s$d


f <- function(numBasis){
  b=0
  for(i in 1:numBasis){
    b = b + d[i]*as.matrix(u[,i])%*%t(as.matrix(v[,i]))
  }  
  b = apply(b, 1, std)
  plot(imagedata(t(b)))  
}

f(10)

f(50)

f(100)

The original picture:

• In this example, eigenbasis is decided by data which is different to other user-defined basis such as cubic spline basis, B-spline basis and wavelet basis.
• We can also use user defined basis so that we only have to store the coefficients in our hard disc rather than the entire image.

Example: image reconstruction (B-spline basis)

library(fda)

g <- function(numBasis){
  bsp = create.bspline.basis(rangeval=c(1,658), norder=3, nbasis=numBasis)
  B1 = getbasismatrix(1:658, bsp, returnMatrix=TRUE)
  bsp = create.bspline.basis(rangeval=c(1,718), norder=3, nbasis=numBasis)
  B2 = getbasismatrix(1:718, bsp, returnMatrix=TRUE)
  D = solve(t(B1)%*%B1)%*%t(B1)%*%y
  A = D%*%B2%*%solve(t(B2)%*%B2)

  est = round(B1%*%A%*%t(B2))
  # est = replace(est, est<0, list(0))
  dim(est) = c(658, 718)
  est.img = imagedata(apply(est, 1, std), type="grey")    
  plot(t(est.img))
}

g(10)

g(50)

g(100)

g(300)

Local polynomial approach

Suppose our goal is to estimate the true function $m(x)$, then, consider local estimation on each data point $x_j$

$$\begin{align*} m(x) \approx& m(x_j) + m'(x_j)(x-x_j) + \cdots + \frac{m^{(p)}}{p!}(x-x_j)^p \\ \approx& \beta_0 + \beta_1(x-x_j) \end{align*}$$ - since we preserve the first two terms of Taylor expansion, $x$ cannot be too far away from $x_j$.

• Now, the parameter of interest is $\beta_0$, we can minimizing the following formulation

$$\min_{\beta_0,\beta_1} \sum_i[y_i-\beta_0-\beta_1(x_i-x_j)]^2$$

• In this formulation, for estimating $m(x_i)$, we still need information over all $x$. To make the Taylor approximate well, we introduce a weight (Kernel function) to emphasize the local information.

$$\min_{\beta_0,\beta_1} \sum_i[y_i-\beta_0-\beta_1(x_i-x)]^2 K_h(x_i-x)$$

• Note that the kernel is a function of bandwidth $h$, it should be prespecified.
• The bandwidth decide how much of the neighboring should be included. It is easily to see that the role of bandwidth is similar to the roughness of penalty in the spline approach.

Local polynomial approach

$$\begin{align*} X=\left(\begin{array}{cc} 1 & x_1-x \\ \vdots & \vdots \\ 1 & x_n-x \end{array}\right), \beta=\left(\begin{array}{c} \beta_0 \\ \beta_1 \end{array}\right), W=\left(\begin{array}{ccc} K_h(x_1-x) & & \\ & \ddots & \\ & & K_h(x_n-x) \end{array}\right) \end{align*}$$ Then, the problem can be reduced to weighted least square $$\min_{\beta}(y-X\beta)^TW(y-X\beta)$$ $$\hat{\beta} = (X^TWX)^{-1}X^TWy$$

• Once we obtain each $\hat{\beta}_0(x_j)$, then connect these points then we got the estimate of function $m(x)$
• Choice of kernel (is not that important compared to the choice of bandwidth)

(ref)

• Choice of bandwidth
  • It is a long story … in simple case, just choose it by your eyes…
• Something we should know
  • compared to spline approach, the development of local theoretical property in local polynomial is much more easier.
  • Something may goes wrong in the boundary.
  • We cannot use smoothing function to do extrapolation.

Curse of dimensionality

• When $p=1$, for axis $x\in[1, 100]$ with integer grid, we can fill this axis with $100$ points.
• When $p=2$, for axis $x, y\in[1, 100]$ with integer grid, we need $100^2$ points to maintain the same density.

• When $p=3$, for axis $x, y, z\in[1, 100]$ with integer grid, we need $100^3$ points to maintain the same density.

• In high-dimensional problem, can you really catch enough neighborhood ?